// Get stopping positions of muons inside the target
#include "common.h"
#include "../rootlibs/TTrackerHit.h"
#include "../rootlibs/TEvent.h"

void GetOrgMomentum(int runNo, TH1F *h);
void ProcessRun(int runNo, TH1F *h)
{
	printf("Run %d:\n",runNo);
	char filename[20];
	sprintf(filename,"run%05d.root",runNo);
	TFile *file = TFile::Open(filename);
	TTree *t = (TTree*)file->Get("t");
	TEvent* evt = new TEvent();
	t->SetBranchAddress("evt",&evt);
	int nevt = t->GetEntries();
	// Real work follows here
	h->SetXTitle("Muon momentum (MeV/c)");
	timeval t1,t2;
	int nstop = 0;
	int nstop0 = 0;
	int nhit_in_target = 0;

	gettimeofday(&t1,NULL);
	for (int i = 0; i < nevt; i++) 
	{
		t->GetEntry(i);
		int nhits = evt->GetNoOfHits();
		for (int j = 0; j < nhits; j++) 
		{
			TTrackerHit aHit = evt->GetHit(j);
			bool is_mu = strcmp(evt->GetHit(j).GetPName(),"mu-") == 0;
			if (is_mu && aHit.GetDetID()==1)
			{
				if(aHit.GetHitType()==1)
				{
					double ke = aHit.GetKE();
					h->Fill(GetMuMomemtum(ke));
				}
				if (aHit.GetHitType()==0)
					nstop++;

				double pos[3];
				aHit.GetPos(pos);
				bool intarget = pos[0]>-25000 && pos[0]<25000 && pos[1]>-25000 
					&& pos[1]<25000;
				if (intarget && aHit.GetHitType()==1)
					nhit_in_target++;
			}	

			if (is_mu && aHit.GetDetID()==0 && aHit.GetHitType()==0)
				nstop0++;
		}
	}
	gettimeofday(&t2,NULL);
	double timelapse;
	timelapse = (t2.tv_sec - t1.tv_sec)*1000000 + (t2.tv_usec - t1.tv_usec);

	std::cout<<"Stop Plastic: "<<double(nstop0)/nevt*100<<std::endl;
	std::cout<<"Stop Al: "<<double(nstop)/nevt*100<<std::endl;
	std::cout<<"hit inside target: "<<double(nhit_in_target)/nevt*100<<std::endl;
	std::cout<<"Time lapse: "<<timelapse<<std::endl;
}

int main (int argc, char *argv[])
{
	TApplication theApp("App",&argc,argv); 
	SetStyle();
	TH1F *h_mu_org;
	h_mu_org = new TH1F("h_mu_org","Muon momentum at target",200,0,40);
	h_mu_org->SetLineColor(28);
	TH1F *h_mu_exit[9];
	for (int i = 0; i < 9; i++) 
	{
		char h_str[20];
		sprintf(h_str, "h%02d",i);
		char h_title[100];
		//int p0 = 30 - i*1;
		sprintf(h_title,"Muon momentum at target, 5 mm plastic scinti.");
		h_mu_exit[i] = new TH1F(h_str,h_title,200,0,33);
		h_mu_exit[i]->SetLineColor(1+i);
	}
	
	//GetOrgMomentum(1,h_mu_org);

	//ProcessRun(1, h_mu_exit[0]); // 0.7 mm
	//ProcessRun(2, h_mu_exit[1]); // 0.8 mm
	//ProcessRun(3, h_mu_exit[2]); // 0.9 mm
	//ProcessRun(4, h_mu_exit[3]); // 1.0 mm
	//ProcessRun(5, h_mu_exit[4]); // 1.0 mm
	//ProcessRun(6, h_mu_exit[5]); // 1.0 mm
	//ProcessRun(7, h_mu_exit[6]); // 1.0 mm
	//ProcessRun(8, h_mu_exit[7]); // 1.0 mm
	//ProcessRun(9, h_mu_exit[8]); // 1.0 mm
	//ProcessRun(10, h_mu_exit[0]); // 0.7 mm
	//ProcessRun(11, h_mu_exit[1]); // 0.8 mm
	//ProcessRun(12, h_mu_exit[2]); // 0.9 mm
	//ProcessRun(13, h_mu_exit[3]); // 1.0 mm
	//ProcessRun(14, h_mu_exit[4]); // 1.0 mm
	//ProcessRun(15, h_mu_exit[5]); // 1.0 mm
	for (int i = 0; i < 6; i++) 
	{
		ProcessRun(10+i,h_mu_exit[i]);
	}

	char canvas_str[50];
	sprintf(canvas_str,"Incoming muon dist");
	TCanvas *c1 = new TCanvas(canvas_str,"All",800,600);
	c1->SetLogy();
	//h_mu_org->Draw();
	TLegend *leg = new TLegend(0.15,0.6,0.31,0.89);
	leg->SetFillStyle(0);
	leg->AddEntry(h_mu_org,"Original muon");
	h_mu_exit[0]->Draw();
	for (int i = 1; i < 6; i++) 
	{
		h_mu_exit[i]->Draw("same");
		char leg_entry[100];
		int p0 = 30 - i*1;
		sprintf(leg_entry,"%d MeV/c",p0);
		leg->AddEntry(h_mu_exit[i],leg_entry);
	}
	leg->Draw();


	//TCanvas *c2 = new TCanvas("c2",canvas_str,900,1200);
	//c2->Divide(2,3);
	//for (int i = 0; i < 6; i++) 
	//{
		//c2->cd(i+1);
		//h_mu_exit[i]->Draw();
	//}

	char pdffile[50];
	sprintf(pdffile,"pdf/Mu_momentum_all.eps");
	c1->Print(pdffile);
	//sprintf(pdffile,"pdf/Mu_momentum_indi.eps");
	//c2->Print(pdffile);

	std::cout<<"Done\n";


	// To keep the ROOT windows open after execution
	theApp.Run(); 
	return 0;
}

void GetOrgMomentum(int runNo, TH1F *h)
{
	printf("Run %d:\n",runNo);
	char filename[20];
	sprintf(filename,"run%05d.root",runNo);
	TFile *file = TFile::Open(filename);
	TTree *t = (TTree*)file->Get("t");
	TEvent* evt = new TEvent();
	t->SetBranchAddress("evt",&evt);
	int nevt = t->GetEntries();
	h->SetXTitle("Muon momentum (MeV/c)");

	for (int i = 0; i < nevt; i++) 
	{
		t->GetEntry(i);
		int nhits = evt->GetNoOfHits();
		for (int j = 0; j < nhits; j++) 
		{
			TTrackerHit aHit = evt->GetHit(j);
			bool is_mu = strcmp(evt->GetHit(j).GetPName(),"mu-") == 0;
			if (is_mu && aHit.GetDetID()==0 && aHit.GetHitType()==1)
			{
				double ke = aHit.GetKE();
				h->Fill(GetMuMomemtum(ke));
			}
		}
	}
}
